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Abstract 

We describe in detail the method for Particle-in cell/Monte-Carlo sim- 
ulation of electron cyclotron resonance (ECR) discharges. In the simula- 
tion, electric and magnetic fields are obtained by solving Maxwell equa- 
tions, and electrons and ions are accelerated by solving equations of mo- 
tion. We consider two different cases: (i) propagation of electromagnetic 
wave in the presence of a constant external magnetic field; (ii) propagation 
of electromagnetic wave in the presence of a linearly decreasing magnetic 
field which corresponds to a realistic ECR discharge. The simulation re- 
sults indicate that at the resonance layer, the electrons are heated by the 
electromagnetic wave, and the incoming wave amplitude is pronouncedly 
damped, with the wave hardly propagating through the ECR layer. 
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In recent years, there has been increasing interest in high density plasmas at 
low gas pressures for semiconductor wafer processing |Jj|. Unlike Reactive Ion 
Etching (RIE) discharges in which the electrons mean free path is of order of a 
few centimeters, in Electron Cyclotron Resonance (ECR) discharges, the elec- 
trons are confined by the external magnetic field and their mean free path is of 
order of less than a millimeter. Thus ECR discharges are capable of generating 
high density plasmas at low gas pressures and low temperatures. Because of 
these advantages, numerous experiments have been performed to study these 
discharges. However, the fundamental understanding of these discharges is not 
yet satisfactory. The main goal of this work is to describe the method for sim- 
ulating ECR discharges using the self consistent Particle-in- cell/Monte-Carlo 
(PIC/MC) technique. 

The present PIC/MC simulator uses particle-in-cell (PIC) scheme for charge- 
assignment- force- interpolation Q , and Monte Carlo technique for collisions and 
scatterings. In the simulation, an electromagnetic wave with frequency luq = 
2.45 GHz enters the system along the z axis. PIC/MC is used to model the 
interaction of the electromagnetic wave with the electrons and ions. 

Very briefly, the PIC /MC algorithm for ECR plasma consists of the following 
five subroutines: 

(J) Interpolate the instantaneous velocities of MC particles representing ions 
and electrons to the grid points using PIC technique. Once the velocity of ions 
and electrons at a grid is obtained, the current density at the same grid can 
easily be calculated. 

(II) Solve the Maxwell's equation on a spatially discretized mesh to obtain 
the electric and magnetic fields. In the simulation, we assume that the electric 
and magnetic fields do not vary along the x and y axes, i.e., we only consider 
E x (z,t), E y (z,t), E z (z,t), B x (z,t), B y (z,t), B z (z,t). This assumption is jus- 
tifiable since electrons are confined to a radius of less than a millimeter by the 
external magnetic field. For example, for an electron temperature of 2 eV, which 
corresponds to the electron velocity of 10 8 cm/sec, and for a magnetic field of 
B = 10~ 2 T, the electron's radius is less than 1 millimeter (here we have used 

the formula r = "^M here 6 is the angle between the momentary electron 
eBsmcp' T to J 

velocity vector and the magnetic field vector, and |v| is the absolute value of 
the electron velocity) . 

Assuming the electric and magnetic fields do not have any variations along 
the x and y directions, the Maxwell's equations can be be written as 

SE * 5B v at n^ 



SB y = SE X 
5t C Sz 



(2) 



SE y 5B X 



st C —-^ ® 

S_B x _S_Ey 
St ~ Sz' W 

f = -4^, (5) 
In the above equations, E's, B's and J's represent the electric field and magnetic 
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field and current density. The boundary conditions are obtained by assuming 
that the circularly polarized waves enter the plasma from the left and that the 
transverse waves leave the system without being reflected. Thus the horizontal 
and vertical components of the electric and magnetic fields at z — and z = L 
are given by 

E x (0) = E cos cat, B x (0) = B sin u>t, E y (0) = E smuit : B y (0) = —B cosut, 
E x (L) = 0, B x (L) = 0, E y (L) = 0, B y (L)=0, (7) 

where w is the frequency of the source, and Eq and Bq are the electric and 
magnetic field of the source, i.e., the electric and magnetic field at z = 0. The 
boundary condition of the longitudinal electric field is obtained by solving Guass 
Equations, i.e., 

E z (0) = 47rtr(0,t), E Z (L) = Ana{L,t), (8) 

where a(0,t) and a(L,t) are the surface charge density at the left and right 
boundaries which vary with time. The surface charge density can be obtained 
from 



<7(0,t)= / (J»(0) - J e (Q))dt', 



a(L,t) = / (J t (L)-J e (L))dt'. (9) 



o 

The longitudinal boundary conditions imply that the total charge of the system 
including the boundaries is zero. 

The Maxwell's equation for the transverse wave is obtained by adding and 
subtracting Eqs. 1 and 2 (or Eqs. 3 and 4), i.e., 

S 5 , 

c— ) F x - y = -AttJ x ., 



St 6 
S__ _6_ 

St C 5z 



G x , y = -4ttJ x;v , (10) 



where F x — E x + B y , G x — E x — B y , F y = E y + B x , and G y = E y — B x . Note 
that the left-hand side of Eq. (10) can be considered as the total derivative 
along the vacuum line z — ct. Thus if we assume Az — cAi, Eq. (10) may be 
discritized as S 

F*»<fi + + = ^ Jx . y (t + ^z + c*) , (11) 

or 

F^l+l = K;y,j - ^ [(J*;vJ+l + J*; V j) /2] At. (12) 

Thus the summation of the x component of the electric field and the y compo- 
nent of the magnetic field at grid point j + 1 and at the time n + 1 depends on 
the summation of the x component of the electric field and the y component of 
the magnetic field at grid point j and at the time n and on the average current 
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density between grids j and j + 1 at time t + At/2. Using the above technique, 
one can obtain E x {z, t), E y (z, t), B x (z, t), and B y (z, t). 

(III) Interpolate the electric field and the magnetic field from the grid points 
to the location of particles. Once the electric and magnetic fields at the location 
of particles are known, equations of motion can be solved. 

(IV) Integrate the equations of motion under the local and instantaneous 
electric and magnetic fields. To move the particles, we have to solve Lorentz 
equation. 



v n+1 = v n + ±At 



E n + \ (v n+1 + v n ) x B r 



m Az y ' 



where g is the magnetic moment and B ex t is the external magnetic field. Note 
that the simulation uses Leap-Frog technique, and hence time n refers to t—At/2 
and time n + 1 refers to t + At/2. 

Since v n+1 appears on both sides of the above equation, one has to proceed 
very carefully. To obtain the velocity at time n + 1 from the velocity at time n, 
we use Boris's technique B, which is based on the following three steps: 
(i) First we define velocity v~ as 

n ,qEAt g AB^t At 
v = v n H t — 14 

m 2 m Az 2 y ' 

Next we define velocity v' which is related to the velocity v~ by the following 
relation 

v' = v~ + v~ x r, (15) 
qB At 

where the function r is defined as r = . 

m 2 

(Hi) Finally we define velocity v+ 

v + = v~ + v' x s (16) 
where the function s is defined as 

2r 

S = IT1VP ,17) 

Boris |0| has shown that the velocity at time t + At can be obtained from the 
following equation 

vn+1 = v+ q _EAt_ L AB^At 

m 2 to Az 2 y J 

The equations of motions for electrons and ions at time t is numerically inte- 
grated to obtain the position of the electrons and ions at time step t + At. 

z(t + At) = z(t) + v n+1 At (19) 

(V) Use random numbers (Monte Carlos technique) and collision cross sections 
to account for scattering and ionizations. The total electron-neutral scattering 
cross section a tota i(v) is <J to tai(v) = ^ t ° tal , where v is the electron velocity 
and Ktotal = 2 x 10~ 8 cm 3 / s is the rate constant. Ionizing collisions occur if 
the electron energy is larger than a specific value (for example, for Argon the 
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threshold energy is 15 eV). An ionizing collision is modeled by loading a new 
electron and ion at the position of the ionizing electron. The kinetic energy 
after ionizing collision is partitioned between the two electrons. Ion-ion charge 
exchange and ion-ion elastic collisions are also included in the simulator. 

In the simulator, the charged particles move under the influence of the self- 
consistent electric and magnetic fields and suffer collisions with neutral particles. 
The neutral gas density is chosen to be 2 x 10~ 14 cm~ 3 . The size of the discharge 
is 24 cm and the number of grids is 2667. Microwave with an amplitude of 0.16 
Gauss and at a frequency of 2.45 GHz enters the system from the left along the 
z axis and propagates through the discharge. 

Figure 1 (a) [Fig. 1 (b)] shows the variations of E X1 E y [B x , B y ] within the 
discharge in the absence of particles. The incoming wave propagates through 
the discharge without attenuation. These figures clearly demonstrate that the 
subroutine solving Maxwell equation is working properly. 

We now consider the coupling of the electro-magnetic wave to the electrons 
and ions. First, we briefly describe the fundamental principles of the ECR 
discharges. We consider an external magnetic field along the z axis with a 
magnitude of B ext (for simplicity we assume that the external electric field 
is zero). An electron rotates around the magnetic field with a frequency of 

lu c = ^jjfj^ (w c is 2.45 GHz when the external magnetic field is B ext = 875 
Gauss) . We now assume that an electromagnetic wave with frequency ojo enters 
this system. If uio is much smaller or much larger than u c , the electron is not 
heated by the incoming wave. However, when luo = lo c , resonance condition 
is attained and the wave energy is absorbed, leading to strong acceleration of 
electrons. 

In the simulation, we assume an electromagnetic wave with frequency ujq = 
2.45 GHz enters the system from the left. First we consider two types of external 
magnetic fields: (1) B ext — 1875 G which corresponds to an electron cyclotron 

frequency of u; c = 3%* = 5.24GHz, (2) B ext = 875 G which corresponds to 
an electron cyclotron frequency of 2A5GHz. Figure 2 (a) [Figure 2 (b)] shows 
the variation of the horizontal component of the electric [magnetic] field within 
the discharge. When B ext = 875 G, which corresponds to electron cyclotron fre- 
quency of 2.45 GHz and which is equal to the frequency of the electromagnetic 
wave, resonance occurs and the electric and magnetic fields of the incoming elec- 
tromagnetic wave are pronouncedly damped, with the wave hardly propagating 
within the discharge. In contrast, when B ext — 1875 G, which corresponds to 
an electron cyclotron frequency of 5.25 GHz and which is much larger than the 
frequency of the incoming wave, the electric and magnetic fields of the electro- 
magnetic wave propagate through the discharge without attenuation (note that 
we use CGS system where both electric and magnetic fields are measured in 
Gauss). 

Figure 3 shows the horizontal component of electron velocity within the dis- 
charge. When B ext = 875 G, resonance occurs and the electrons are heated 
by the incoming electromagnetic wave. Thus the horizontal component of the 
electron velocity increases sharply right at the boundary. However, When 
B ext — 1875 G, the electrons do not absorb much energy from the incoming 
wave. The electron velocity is therefore small at the boundary and does not 
change rapidly within the discharge. 

In the previous examples, we assumed that the external magnetic field is 
constant within the discharge. In an actual ECR discharge, B ext drops along z 
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axis. To model a realistic ECR discharge ||, we simulated a system where the 
external magnetic field decreases linearly along the z axis so that at the center of 
the discharge, the external magnetic field is 875 G and hence resonance occurs. 

Figs. 4 (a) and 4 (b) show the trajectories of electrons at the center of 
discharge when the external magnetic fields of 1875 G and 875 G. Although in 
the simulation, we trace the electrons only along the z direction, Figs. 4 (a) 
and 4 (b) are obtained by integrating V x and V y over time. Note that when 
B e xt = 875 G, i.e., when uj c = ujq, resonance condition is attained and electrons 
spiral around the external magnetic field. In contrast, when B ext = 1875 G, i.e., 
when lu c >> luq, resonance does not occur and electrons do not spiral around 
the magnetic field. 

Figures 5(a) and 5(6) show the horizontal and vertical components of the 
electric field within the discharge. The incoming electromagnetic wave is en- 
tirely absorbed by the electrons near the resonance layer, with both E x and E y 
dropping rapidly near the ECR layer. The incoming wave hardly propagates 
beyond the resonance layer. 

Figures 6(a) and 6(6) show the horizontal and vertical components of the 
electron velocity within the discharge after 65 cycles. At the resonance layer, 
the velocity of electrons jumps significantly, indicating that the electrons ab- 
sorb considerable energy from the incoming electromagnetic wave. However at 
other positions within the discharge, the velocity of electrons does not change 
rapidly, indicating that the electrons absorb very little energy. Of course after 
many cycles, the horizontal and vertical components of electron velocity be- 
come isotropic as the the electrons localized near the resonance layer create new 
electrons by ionization as well as diffusing toward the boundaries. To clearly 
demonstrate the heating of the electrons at the resonance layer, we present here 
the simulation results after only 65 cycles. 

In summary, PIC/MC technique has been used to investigate the funda- 
mental properties of an electropositive ECR discharge. The simulation results 
indicate that at resonance layer, i.e., at oj c — loq, the incoming electromagnetic 
wave is pronouncedly damped, leading to dramatic acceleration of electrons. 
The simulation results also show that at resonance layer, electrons spiral around 
the external magnetic field. PIC/MC technique shows great promise for sim- 
ulating more complex (for example electronegative) discharges in two or three 
dimensions. 
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Figure Captions 

IMPORTANT NOTE: 

Figs. 1 (a) and 1 (b) are on the same page. 
Figs. 2 (a) and 2 (b) are on the same page. 
Figure 3 is on one page. 
Figs. 4 (a) and 4 (b) are on the same page. 
Figs. 5 (a) and 5 (b) are on the same page. 
Figs. 6 (a) and 6 (b) are on the same page. 

Fig. 1 (a) Profile of the horizontal (solid line) and vertical (dashed line) com- 
ponents of the electric field within the discharge in the absence of particles. 
Fig. 1 (b) Profile of the horizontal (solid line) and vertical (dashed line) com- 
ponents of the magnetic field within the discharge in the absence of particles. 
Fig. 2 (a) Profile of the horizontal component of the electric field within the 
discharge when the external magnetic field is at 875 G (solid line) and at 1875 
G (dashed line). 

Fig. 2 (a) Profile of the horizontal component of the magnetic field within the 
discharge when the external magnetic field is at 875 G (solid line) and at 1875 
G (dashed line). 

Fig. 3 Profile of the horizontal component of electron velocity within the dis- 
charge when the external magnetic field is at 875 G (solid line) and at 1875 G 
(dashed line). 

Fig. 4 (a) Trajectories of electrons when the external magnetic field is at 1875 
Gauss. 

Fig. 4 (a) Trajectories of electrons when the external magnetic field is at 875 
Gauss. 

Fig. 5 (a) Profile of the horizontal component of electric field in an ECR dis- 
charge. Note that at the center of the discharge, the external magnetic field is 
875 G. 

Fig. 5 (b) Profile of the vertical component of electric field in an ECR discharge. 
Fig. 6 (a) Profile of the horizontal component of electron velocity in an ECR 
discharge after 65 cycles. 

Fig. 6 (b) Profile of the vertical component of electron velocity in an ECR 
discharge after 65 cycles. 
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Horizontal component of the magnetic field (Gauss) Horizontal component of the electric field (Gauss) 

-0.15 -0.05 0.05 0.15 -0.1 0.0 0.1 



Horizontal electron velocity (cm/sec) 



-6 -5 -4 -3 -2 -1 

e,8 e,8 e,8 e,8 e,8 e,8 e,0 
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Vertical electric field (Gauss) Horizontal electric field (Gauss) 

-0.10 0.0 0.05 -0.10 0.0 0.05 




